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ABSTRACT 

> ■ 

\Q . This paper explores orbits in extended mass distributions and develops an 

analytic approximation scheme based on epicycloids (spirograph patterns). We 
focus on the Hernquist potential if) = 1/(1 + £), which provides a good model 
for many astrophysical systems, including elliptical galaxies (with an i? 1 / 4 law), 
dark matter halos (where N-body simulations indicate a nearly universal density 
^! profile), and young embedded star clusters (with gas density p ~ £ _1 ). For a given 

potential, one can readily calculate orbital solutions as a function of energy and 
angular momentum using numerical methods. In contrast, this paper presents 
a number of analytic results for the Hernquist potential and proves a series of 
general constraints showing that orbits have similar properties for any extended 
mass distribution (including, e.g., the NFW profile). We discuss circular orbits, 
radial orbits, zero energy orbits, different definitions of eccentricity, analogs of 
Kepler's law, the definition of orbital elements, and the relation of these orbits to 
spirograph patterns (epicycloids). Over a large portion of parameter space, the 
orbits can be adequately described (with accuracy better than 10%) using the 
parametric equations of epicycloids, thereby providing an analytic description of 
the orbits. As an application of this formal development, we find a solution for 
the orbit of the Large Magellanic Cloud in the potential of our Galaxy. 



Subject headings: stellar dynamics, celestial mechanics — methods: analytical 
— galaxies: kinematics and dynamics, halos — stars: formation 
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1. INTRODUCTION 

Orbits are a fundamental component of dynamical astrophysical systems, including 
dark matter halos, galaxies, star clusters, and solar systems. For a single point mass, the 
source term for a Keplerian potential, orbits are well-known conic sections and have been 
the cornerstone of solar system dynamics for centuries (starting with Newton 1687; for a 
modern treatment, see Murray & Dermott 1999). Many astrophysical systems have extended 
mass distributions, which give rise to more complex potentials and more complicated orbits 
(Binney & Tremaine 1987). Although an enormous amount of work has been done on 
understanding the dynamics of such extended mass distributions, the orbits themselves are 
not as well studied. In this contribution, we calculate the orbits for a class of extended mass 
distributions that provide good working models for dark matter halos, elliptical galaxies, 
and embedded young star clusters, albeit in rather different regimes of parameter space (see 
also van den Bosch et al. 1999 for a consideration of orbits in an isothermal potential). 
Along with our improved characterization of the orbits, we develop an effective analytic 
approximation scheme that describes the orbits in terms of epicycloids (spirograph patterns) 
and thereby obtain a set of parametric equations for the orbits. In addition to making 
dynamical calculations easier, this analytic approximation scheme adds to our conceptual 
understanding of the orbits — in much the same way that knowing Keplerian orbits are 
conic sections provides increased perception. 

In this paper, we focus on the Hernquist potential, which arises from a density profile 
of the form 

where the dimensionless radius £ = r/r s and r s is the length scale of the potential. This 
density profile has a corresponding analytic form for its distribution function (Hernquist 
1990). One goal of this paper is to characterize the orbits for bodies traveling within this 
extended density distribution. The motivation for this effort is that particular subclasses 
of these orbits describe the motion of bodies in many astrophysical systems. The original 
application of this potential (Hernquist 1990) was to approximate the i? 1//4 law for elliptical 
galaxies (de Vaucouleurs 1948). In addition, density profiles of this form appear in dark 
matter halos of galaxies and galaxy clusters, and in the molecular cloud cores that form 
stellar groups and clusters. Both of these latter applications require some explanation: 

Over the past decade (starting with Navarro, Frenk, & White 1997; hereafter NFW), nu- 
merical simulations of cosmic structure formation have shown that the density distributions 
of dark matter halos approach a nearly universal form (see also Crone et al. 1996, Moore et 
al. 1998, Bullock et al. 2001). Although the original form for this universal profile (NFW) is 
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not as steep as that of equation (1) at large radii, subsequent work indicates that the Hern- 
quist density profile provides a good description for the asymptotic structure of dark matter 
halos (Busha et al. 2003, 2004). One advantage of using equation (1) is that the total mass 
of the system is finite (whereas the integral of the original NFW density profile diverges). 
One long term goal of structure formation studies is to understand the dynamics of these 
halo systems, including both the origin of the universal profile and its consequences. Toward 
that end, we must describe the orbital motions of test particles moving in the gravitational 
potential created by the universal density profile (where the "test particles" can range from 
individual dark matter particles with mass mw ~ 100 GeV to small satellite galaxies with 
mass M sat ~ 10 9 M ). 

Molecular cloud cores that form stellar clusters provide another application. They are 
inferred to have equations of state softer than isothermal (e.g., Larson 1985; Jijina et al. 
1999) and hence density profiles of the approximate form p ~ 1/r, which coincides with the 
Hernquist density profile at small radii. The orbits found for the inner limit of the density 
profile (1) are thus applicable to star forming regions. Most stars form in groups/clusters 
(Lada & Lada 2003, Porras et al. 2003) and subsequent stellar orbital motions determine 
the degree to which the environment affects star and planet formation (e.g., Adams & Myers 
2001; Adams & Laughlin 2001). 

Since the potentials considered here are spherically symmetric, motion is confined to a 
plane (in the absence of perturbations) and the orbits can be calculated to high accuracy 
using numerical methods. This paper takes a complementary approach and determines the 
nature of the orbits using analytic methods as far as possible and develops an analytic 
approximation scheme based on spirograph patterns. Our goal is to understand the orbits at 
the same level that we now understand Keplerian orbits, in spite of the greater mathematical 
complexity. Several subtleties arise that cannot be captured through numerical solutions 
alone: The orbital eccentricity can be defined in several ways, and each definition agrees 
with the Keplerian result in different aspects. The orbits do not close and their turning 
angles satisfy A0 < 7r (for a half orbit) and this angular deficit remains even in the circular 
limit. The maximum angular momentum, the maximum turning angle, and the radius of the 
circular orbit can all be determined analytically as a function of energy. These results allow 
for many properties of these orbits to be found analytically and add to our understanding 
of how the orbits depend on the underlying potential (see §2). 

Although a large number of results can be obtained analytically, a complete description 
of the orbital shapes requires numerical evaluation. Nonetheless, the orbit shapes are similar 
to the familiar form of epicycloids, often known as spirograph patterns (and not to be con- 
fused with epicycles - see below). This paper takes this analogy one step further (§3) and 
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shows that for a large fraction of the relevant parameter space, the orbits in the Hernquist po- 
tential can be modeled with high accuracy (better than 10%) using the parametric equations 
of epicycloids. The resulting spirographic approximation provides an analytic description of 
the orbital shape. This analytic description, in turn, allows for other orbital quantities to 
be calculated analytically. For example, we show the effects of geometry on orbital decay by 
integrating Aj/j over an orbit (§3.3) and we use the spirographic approximation to model 
the orbit of the Large Magellanic Cloud (LMC) in the potential of the Galaxy (§4). 

In addition to the applications presented in this paper, the analytic approximations 
developed here can be used in a wide variety of other contexts. For example, an analytic 
description of the orbits can be useful in galactic stability theory (e.g., Evans & Read 1998), 
where the perturbed distribution function is given as an integral over the equilibrium orbit. 
Another application of analytic orbits is to find steady-state distribution functions, which 
depend on the isolating integrals of motion (or the constants in the spirographic approxi- 
mation developed here). This analytic description of orbits will also be useful in assessing 
the effects of cluster environments on the star formation process: In this context, the orbits 
of newly formed stars determine how much radiation impinges on their circumstellar disks 
and their probability of experiencing scattering interactions; this physical processes, in turn, 
determine how the cluster environment can compromise planet formation. 

Although extended mass distributions do not always have the form of equation (1), we 
argue that the orbits found here provide good baseline models for orbital behavior in more 
general potentials. In particular, generic orbits in extended mass distributions will be more 
like those of the Hernquist model than the conic sections resulting from a Keplerian potential. 
Further, such generic orbits can be described through the spirographic approximation devel- 
oped in this paper. Orbits in general extended mass distributions share a number of traits: 
For mass distributions of physical interest, the topology of the manifold of orbital solutions is 
simple (the same as the Si x Si structure of the Kepler problem; see Smale 1970ab), although 
the geometry of the orbits is more complicated. One defining characteristic of these orbits - 
turning angles for a half orbit bounded by tt/2 < A6 < tt - holds for general extended mass 
distributions (Contopoulos 1954). This paper proves additional results concerning generic 
orbits in extended mass distributions (§5). We show that the eccentricity can be defined in 
multiple ways, with the definitions coincident only for a Keplerian potential. We also show 
that orbital shapes in generalized extended mass distributions are confined to be relatively 
close to those of epicycloids (spirograph patterns). We validate this claim by finding a set of 
bounds on orbital shapes that apply to all physical potentials, and by explicit construction 
of orbits - and comparison with the spirographic approximation - for a collection of well 
known potentials. 
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2. ORBIT SOLUTIONS 

In this section we calculate orbits for bodies moving in the potential of the Hernquist 
profile (with the density distribution of eq. [1]). We define the total depth of the gravitational 
potential well = 2TrGp r 2 an d the total mass = 27rrfp , so that the corresponding 
mass profile, force profile, and gravitational potential take the forms 

M(0 = M -(iT^' F(0= (TT^' m = ITr (2) 

All quantities are taken to be positive (with the proper signs inserted as necessary). 



1 -1/2 

2(E - V)r 2 /f - 1 (3) 



2.1. General Orbits 

The basic formulation for orbits in spherical mass distributions is well known (Binney & 
Tremaine 1987); this section defines notation and presents analytic results for orbits in the 
Hernquist potential. An orbit in any potential V(r) can be determined from the differential 
equation 

d9 _ 1 
dr r 

where E is the energy and j is the specific angular momentum. This paper focuses on bound 
orbits with negative energy. We define dimensionless variables for the energy (e) and angular 
momentum (q) according to 

e=\E\/* and q = f/2V r 2 s . (4) 

The differential form of the orbit equation can then be written t^dO/dt^ = [q/ f{i)] 1 ^ 2 , so that 
much of the dynamics is determined by the properties of the cubic equation 

/(0 = -^ 2 + I ^|-9- (5) 

For a given energy e and angular momentum parameter q, the orbit exists within the range 
of £ for which /(£) is positive. The two positive zeroes of /(£) thus correspond to the radial 
turning points of the orbit and are denoted here as £i and £ 2 . The third zero of the cubic 
is always negative for the physical range of parameter space and we denote this root as 
a = — £ > 0. Finding the roots of the equation (5) is straightforward, but cumbersome 
(e.g., Beyer 1980). However, the inverse relations can be found simply: 

6 + 6 + 66 (66) 2 



' = (6 + 6) (i + 6 + 6 + 66) ' Q (6 + 6) (l + 6 + 6 + 66) ' 
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and a = -£ = ~, j^r ~, ~, ■ (6) 

One quantity of interest for a given orbit is the angle A9 through which the orbit turns 
over one half cycle, e.g., as the orbit moves from £1 to £ 2 . This angle is given by 

M = Vq rf[/(0r V2 - (7) 

A related quantity of interest is the orbital half period, denoted here as T\/2- The natural 
time scale of this potential is r = r s (2^ )' 1 ^ 2 . The dimensionless time scale Ti/ 2 = ti/ 2 /t 
can be written in integral form 

ri/2= r^e[/(0]" 1/2 - (8) 

For a given dimensionless energy e, orbits can have a maximum specific angular mo- 
mentum and hence a maximum value of q. This maximum value takes the form 

1 (1 + Vl+We - 4e) 3 

8e (l + v /T+8^) 2 ' [) 

In the limit of low energy e — > 0, we recover the Keplerian result 4eg ma:r = 1. In the opposite 
limit of e — > 1, we write e — 1 — 5 and find the limiting form q max = 4<5 3 /27. 

For a given energy e, the minimum of the cubic function /(£) is independent of the 
angular momentum variable q. The location of this minimum, denoted here by where 
df /<%{£*) = 0, is given by 

In the limit of low energy e — > 0, we recover the Keplerian form ^ = l/2e. As a result, we 
can identify as playing the same role that the semi-major axis plays for bound orbits in 
Keplerian potentials (although these orbits are not ellipses). Notice also that is the radius 
of the circular orbit for a given energy e. 

For given values of energy and angular momentum (e, q), the orbit takes a specified 
form. One such orbit is shown in Figure 1 for the case e = 0.10 and q = q ma x/2- For each 
half oscillation in the radial direction, as £ ranges from £ 2 down to £i, the orbit sweeps 
through an angle A9, which lies in the range tc/2 < A9 < n. In the Keplerian limit (e — > 0), 
A# — > 7r. In the limit of radial orbits (q — > 0), the angle A6 1 — > n/2. Between these limits, 
the angular displacement takes on intermediate values that depend on energy and angular 
momentum. These results are illustrated in Figure 2, which shows A# as a function of 
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q/q m ax for a collection of energy values. One curious result is that the angle A9 does not 
approach n in the limit of circular orbits (see eq. [16] below). For the set of curves shown 
in Figure 2, specifically for energy and angular momentum in the range 0.2 < e < 1 and 
1CT 6 < q/q ma x < 1, the turning angle A9 can be approximated with a fitting function of the 
form 

— = - + [(1 + 86)- 1 / 4 - + ln(g / W) ] 3 - 6 . (11) 
7T 2 LV ; 2 JL 6 In 10 J K J 

This function provides an approximation to the turning angle with an accuracy of about 
2.5% over the previously stated portion of parameter space. 

Figure 3 shows the orbital half period, i.e., the time required for the orbiting body to 
travel from the inner turning point to the outer turning point. In the figure, the half period 
is normalized by the factor e 3//2 , which represents the variation expected in the Keplerian 
limit, and by the factor [cos -1 \fe — \/i(l — e) 1 / 2 ] that gives the correct orbit time in the limit 
of zero angular momentum. The resulting (normalized) half period is nearly independent 
of angular momentum, and nearly the same as that of the radial orbit. The orbiting body 
spends most of its time at large £, where q has little effect, so the result is almost independent 
of q. The half period does, however, depend on the energy e. 

Given an energy and angular momentum (e, q), and hence the radial turning points, one 
can define a generalized "eccentricity" in terms of £i and £ 2 . Unlike the Keplerian case, the 
definition of eccentricity leads to some ambiguities. Two reasonable definitions of eccentricity 
are 

_ £2 - £1 , ~_£2-£i / 1oN 

e= — and e = — . (12) 

At low energies e — » 0, the orbits become Keplerian and both definitions of the eccentricity 
coincide with the traditional one. In general, however, the two definitions do not coincide 
(see §5.3). The first definition of eccentricity e has the advantage in that it includes the 
radius £* of the circular orbit and the turning points represent symmetric departures from 
(analogous to the Keplerian case). In contrast, the second definition e does not correctly 
describe small departures about the circular orbit because £2 + £i 7^ 2^*. On the other hand, 
the second definition allows for the full range of eccentricities < e < 1, whereas the first 
definition has a maximum eccentricity so that < e < e max (for all nonzero energies e > 0). 
This maximum eccentricity occurs for radial orbits (q — > 0), where the outer turning point 
£2 — > (1 — e)/e and the inner turning point £1 — > 0, and is given by 

e - 2(1 ~ C) (13) 

In the Keplerian limit e — > 0, one recovers the expected result that e max — > 1. In the inner 
limit e — > 1, one finds e max — > 3/4. In between, the maximum eccentricity is a continuous 
monotonic function of dimensionless energy e. 
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2.2. Limiting Forms of the Solutions 

For the cases of circular orbits, radial orbits, and zero energy orbits, the solutions reduce 
to limiting forms that can be readily evaluated. For circular orbits, one can solve for the 
relationship between the orbital period r and radius. The resulting "Kepler Law" for circular 
orbits takes the form 

- 2 = ^(l + 2 - (14) 

For a circular orbit at radius £, the energy and angular momentum are given by e = (1 + 
£/2)(l + £)~ 2 an d Q = (1/2)£ 3 (1 +£)~ 2 - A s a re sult, we can write the orbital period in terms 
of the dimensionless energy so that 



r 



2 = jn_ 1 

Gp e 



1 + 46+^1+86 

8? 



" l] ■ (15) 



As defined here, the period r is the time required for the orbit to make a complete circuit, 
i.e., for the angle to trace through 2n radians. However, the angular integral (eq. [7]) does 
not approach n in the limit of circular orbits (as it does for Keplerian orbits). The angular 
displacement is less than it and depends on energy according to 

lim A# = 7r(l + 8e)- 1/4 . (16) 

In the Keplerian limit e — > 0, A0 — > 7r, as expected. For f ^ 0, the formula (16) agrees 
with the limit points of the curves shown in Figure 2. Since circular orbits are stable in this 
potential, small perturbations about the circular orbit lead to oscillations at the epicyclic 
frequency re, which takes the form re 2 = 2irGp (3 + + £) -3 - 

Another class of orbits that can be considered separately are those with zero angular 
momentum, i.e., radial orbits. In this case, the angular coordinate is no longer of interest 
and the orbits are characterized by their energy e. Each energy has a corresponding infall 
time scale r, defined to be the time required to fall from the outer turning point £2 to the 
center. This time scale can be evaluated to obtain the result 

r = (4nGp )- 1/2 e- 3 / 2 [cos" 1 + vWT^l . (17) 

For radial orbits, the position along the orbit is specified by the radial coordinate so the 
angle is no longer necessary. Nonetheless, the angle of the orbit remains defined and the 
half angle of the orbit (eq. 7) reaches a definite limit. As q — > 0, we find A# — > tt/2. 

Finally, we consider zero energy orbits. For Keplerian potentials, such orbits are parabo- 
las. Here, the corresponding orbits are more complicated, but simple solutions can be found 
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for limiting cases. The key parameter for a zero energy orbit is its inner turning point, which 
is given by 



6 = 



q + ^/q 2 + 4q 



(18) 



One limiting form of this solution occurs when £ <C 1 for the entire orbit. This solution 
would also result for the case in which a density profile of the form p ~ 1/r extends out 
to infinity. In this limit, the orbit equation can be solved to obtain sin6> = y/q/£, i.e., the 
orbital path is a straight line defined by y = £ sin 6 = y/q = constant. For orbits that remain 
in the outer halo where £ ^> 1, we can also simplify the expression to obtain 



cos 



^ + 1 x1/2 

' q 



2(g+l) 



(19) 



2.3. Orbital Elements 

A particle in orbit can be described by six phase space variables, which are usually 
taken to be the position vector x = (x, y, z) and the velocity vector v = (v X) v y , v z ). For the 
Kepler problem, one can replace these phase space variables with osculating orbital elements, 
i.e., the elliptical (Keplerian) orbit that has the same instantaneous position and velocity 
vectors (Murray & Dermott 1999). For trajectories in extended mass distributions, the phase 
space variables can be written in terms of the elements of the orbits found here. Since the 
potential is spherically symmetric and angular momentum is conserved, orbital motion is 
confined to a plane. The two angles that define the unit normal of the orbital plane (9 P , (f> p ) 
thus provide the first two orbital elements. The rest of this discussion is confined to this 
orbital plane. Although the assignment of the remaining orbital elements is not unique, we 
present a sensible prescription in what follows. 

Within the plane, the orbital shape is determined by the energy e and angular momen- 
tum, specified here by q. The pair (e, q) thus can be used as the next two orbital elements. 
However, other choices are available: The energy e can be replaced with the effective semi- 
major axis £*, as defined by equation (10). For a given energy and angular momentum, the 
cubic equation (5) defines the turning points £i and £2 of the orbit, so that the generalized 
eccentricity (eq. [12]) is well defined and can be used as the other variable. One can thus 
replace the variables (e, q) with the pair (£*,e), which provides the closest analogy to the 
Kepler problem. Alternately, one can use the turning points themselves, (£1,^2), to specify 
the orbital shape. Recall that the cubic equation that defines the turning points has three 
real roots in the regime of interest; because the problem has only two parameters (e,q), 
specification of the turning points (£1,^2) is sufficient to define the equation, and the third 
root can be written in terms of the turning points (eq. [6]). 



-10- 



The two remaining orbital elements must specify the orientation of the orbit within 
the plane and the position of the particle along its orbit. In the custom of solar system 
dynamics, the term longitude is used to denote an angle that is measured relative to a fixed 
location in inertial space. Here we can specify the orientation of the orbit in terms of the 
longitude w of the inner turning point £1 (the longitude of pericenter). Like the elliptical 
orbits in the Kepler problem, the longitude of pericenter can precess. Unlike the case of 
the Kepler problem, orbits in extended mass distributions are not closed and the longitude 
zu of pericenter changes with every orbit. A full orbit turns through an angle 2A9 < 2ir, 
so the pericenter precesses backwards at the rate of 2(n — A9) per orbit. The longitude of 
pericenter is given by w = 2(n — A9)n + zu , where zu is the starting value (t = 0) and n 
is the number of full orbits that have elapsed. Notice that we can use either the starting 
longitude Wq or the current longitude of pericenter w as the next orbital element. For 
completeness, we can also define a continuous longitude of pericenter through the relation 
w c = 2(n — A6)t/r + wq, where r is the (full) orbital period. Finally, the location of the 
particle along the orbit is given by its angular displacement 9. Again following the custom 
of solar system dynamics, the angular displacement variable can be called the anomaly. The 
angle 9 is thus the true anomaly and can be measured in two different ways. The natural 
definition is to measure 9 from the inner turning point (or the outer turning point) for the 
current longitude of pericenter w. However, one could also measure the angle with respect 
to the starting longitude of pericenter zuo , in which case the anomaly is denoted as 9$. We 
can also define a mean anomaly M = uj(At), where At is the time since pericenter passage 
and the mean motion variable uj = (A0)/ti/ 2 (see eqs. [7-8]). 

A collection of possible orbital elements is summarized in Table 1. Here we assume 
that the orbital plane is specified so that only four phase space variables remain. The 
usual variables [(x, y), (v x , v y )\ can be replaced with various choices for the osculating orbital 
elements listed in Table 1. For each pair of shape parameters, any of the possible pairs for 
the orientation and anomaly are viable choices. The last line of the table lists the variables 
required to describe the orbit as an epicycloid (spirograph pattern), as discussed in the 
following section. 

Table 1: Orbital Elements 



Shape Parameters 


(Orientation, Anomaly) 


(e, ?) 




(£*, e) 


(w,9) 


(6, 6) 


(w c , M) 


(a,/3,7) 


(0, i p ) 
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3. EQUIVALENT SPIROGRAPHIC ORBITS 

The orbits found here for the Hernquist potential look much like epicycloids, more 
generally known as spirograph patterns. In this section, we explore the relationship between 
the actual orbits and the closest equivalent spirographic orbits. We find a quantitative 
measure of the difference between the physical orbits and the spirographic approximation, 
and show that the spirographic treatment is valid over a large portion of parameter space. 
This finding allows for a class of physical orbits to be described using simple parametric - and 
hence analytic - equations. The resulting analytic description facilitates a greater physical 
understanding of the orbital dynamics. For example, the effects of dissipation (friction) on 
the orbits can be modeled analytically using the spirographic approximation (see §3.3). 

The epicycloids discussed here are not the same as epicycles, which are often considered 
in connection with astronomical orbits. In the original Ptolemaic definition, an epicycle is 
the curve generated by a small circle that rotates while it moves along the circumference 
of a larger circle. Epicycloids are curves generated by a smaller circle rotating inside (or 
outside) a larger circle, with a drawing point located a given distance away from the center 
of the smaller circle. In modern astronomy, the epicyclic frequency is the frequency of 
radial oscillation for a body in a circular orbit, when perturbed in the radial direction 
by an infinitesimal amplitude [so that k 2 = r~ 3 <i(r 4 f2 2 ) / 'dr\. As a result, the inclusion of 
epicyclic motion (radial oscillations with frequency k) provides a first order correction to a 
circular orbit in an arbitrary potential [which determines the rotation curve fi(r)]. Many 
previous studies have explored the epicyclic approximation and the inclusion of higher order 
terms (e.g., Shu 1969, Kalnajs 1979, Dehnen 1999). In contrast, this present work explores 
epicycloids as approximations to physical orbits. 

3.1. Parametric Description of the Orbits 

In a fixed plane, here the orbital plane, the parametric equations for an epicycloid (a 
generalized spirograph pattern) can be written in the form 

x(i p ) — (a — (3) cosi p + 7cos[(o; — (3)i p / (3] , 

y(t p ) =-(«-/?) sint p + 7 sin[(a - (3%/ (3} , (20) 

where a is the radius of the larger spirograph circle, (3 is the radius of the smaller circle, 
and 7 is the distance from the center of the smaller circle to the drawing point (e.g., see the 
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website of Little 1997, 1 which includes interactive drawing capability). These length scales 
should be considered dimensionless; they can be converted into physical units by specifying 
a length scale (r s ). Keep in mind that the variable t p is a parameter and does not represent 
physical time (its relationship to time is elucidated below). With this form of the parametric 
equations, the spirograph pattern begins at t p — with its maximum radius with an angle 
6 = arctan(y/a;) = 0; the pattern then turns (counterclockwise) through an angle (1 — (3/a)ir 
as the radius decreases to its minimum value. The parametric variable t p ranges from to 
fin I a over this half orbit. Notice that the starting orientation of the orbit can be changed by 
including a starting phase angle (equivalent to an offset in the parametric time variable). 

The parametric equations (20) allow for two different (equivalent) representations of 
physical orbits (see below). In this treatment, we use the representation where 7 plays 
the role of the semi-major axis. In this case, the transformation between the spirographic 
variables and the physical variables - and the inverse transformation - take the form 

& = 7 - (a - P) , £2 = 7 + (a - 0) , A6 — (a — f3)n/a , 

7 = ^(6 + 6), a = ^b^-£i), /3 = i[7r/(Atf)-i](6-6). (21) 

Since the turning points £1 and £ 2 of the physical problem are related to the energy e and 
angular momentum (specified by q) through equation (5), the above relations (eq. [21]) 
implicitly define the spirograph parameters in terms of the energy and angular momentum. 
The physical problem contains only two variables, which can be taken to be the energy and 
angular momentum (e, q) or the turning points (£1, £ 2 ). The equivalent epicycloid is written 
in terms of three parameters (a, (3, 7), but only two are independent. Any such triple will 
define an epicycloid, but only those related through equations (21) are physically relevant. 
The spirographic parameters are related to the energy and angular momentum of the orbit 
via 

e 7 2 + 27 -(a-l) 2 [7 2 -(tt-D 2 ] 2 m) 

2 7 [(l+ 7 ) 2 -(«-/3) 2 ] q 2 7 [(l + 7) 2 -(«-/9) 2 ]' 1 ) 

Notice that since £1, £2 = 7 ± (ot — f3), the parameter 7 is akin to the semi-major axis of 
the orbit (although this analogy is not exact because the orbits are not ellipses). A natural 
definition of the eccentricity e s for a spirographic orbit has the form 

e s = . (23) 

7 

Since 7 > {a — ($) for physically relevant orbits, the eccentricity e s < 1. 



1 Little, D. P. 1997, http://www.math.dartmouth.edu/dlittle/java/SpiroGraph/ 
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In the Keplerian limit, the turning angle A6 — > ir. Thus, the Kepler limit corresponds 
to the limit in which /3 < a and the turning points are given by £ 12 = 7 ± a. In the limit of 
circular orbits, the above expressions (eqs. [22 and 23]) indicate that (a — f3) — > 0. However, 
in order for the parametric equations (20) to trace a circle, the condition (a — (3)/ (3 — > 1 
(a/ (3 — > 2) must also hold. In order for both of these limits to apply, a,f3 — > 0. The limit of 
radial orbits is realized when the angular momentum vanishes, which requires 7 = (a — (3), 
and A9 — > ir/2, which requires a = 2f3\ as a result, radial orbits are given by (a, (3, 7) = 
(2x, x, x) for any value of x. 

The range of possible epicycloid patterns is much larger than the range of orbital shapes 
that represent physical motion (in the potential of an extended mass distribution). The 
allowed parameter space is depicted in Figure 4 (where the parameter 7 is set equal to unity) . 
For physical orbits, the turning angle for a half orbit falls in the range n/2 < A9 < n; for our 
chosen representation, this constraint restricts the parameter space to the region (3 < a/2. 
Since the inner turning point must be greater than zero, 7 > (a — f3). The allowed region of 
the a — (3 plane is thus bounded from above by the dashed line and from below by the dotted 
line. Above the dashed line where f3 = a/2, the turning angle is too small and the resulting 
spirograph patterns are more open than those realized in physical potentials. Below the 
dotted line, the inner turning point is formally negative, so that the pattern turns "in front 
of" the center, rather than orbiting around it. The resulting pattern is tighter than those 
realized in physical orbits. The inset diagrams show the basic orbital shapes for the different 
regions of parameter space. Notice that for this formulation, 7 > a, j3, and the drawing 
point of these epicycloids extends beyond the smaller circle (so they cannot be drawn with 
conventional spirographic wheels). 

The region of the a — (3 plane above the dashed line and below the dotted line (the 
upper right portion of Figure 4) corresponds to the second representation of the orbits. In 
this region, the complementary epicycloid pattern has the same shape, turning points, and 
turning angle as the "original" pattern under the transformation 

y = (a-/3), (a'-/3')=7, a' / (3' = (1 - (3/a)-\ t' p = -{a/ (3 - 1% , (24) 

where the primes denote the complementary variables. The sign transformation of the para- 
metric time variable is necessary to make the spirographic pattern turn in the counterclock- 
wise sense; the factor {a/ (3 — 1) is not necessary but it keeps parametric time passing at 
the same rate. For completeness, we note that epicycloids with (3 < are also allowed and 
correspond to the smaller circle rotating on the outside of the larger circle (although such 
solutions are not considered here). 

The parametric equations (20) completely specify the geometry of the orbit (within 
the spirographic approximation). The velocity is given by the derivatives of the parametric 



-14- 



equations with respect to physical time, i.e., 

v x = -{(a - (5) sintp + 7 (a//3 - 1) sin [(a - P)i P /P]}^, 



v y = | -(a - (5) cosi p + 7 (a//3 - 1) cos[(a - f3)i p /f3]^^ . 



(25) 



The direction of the velocity vector is thus determined (independent of the derivative di p j dt) 
for any point along the orbit. Since the radial displacement is given by £ = (x 2 + y 2 ) 1 ^ 2 , the 



magnitude of the velocity is determined from conservation of energy, which can be written 
in the dimensionless form 

(26) 



2 2,2 

V =V X +V y 



1 + e 



where v 2 (physical) = v 2 (2^ ). We can complete the description by using this result to find 
the variation of parametric time t p with physical time t, i.e., 



dt p 
~dt 



r 1 


1/2 - 







! 2 {a/f3 - l)a/P +{a- (3) 2 a/f3 - {a/f3 - 1)£ 2 



-1/2 



(27) 



This expression gives dt p /dt in dimensionless form. This formalism can be converted into 
physical units, where di p /dt has units of inverse time, using the fiducial time scale defined 

tyfg = (r./2tf ). 

The full orbits in this problem sweep through the angular coordinate faster when the 
orbiting body is near the inner turning point (and slower near the outer turning point). 
In this approximation scheme, both the parametric expressions (eqs. [20] and [25]) and the 
conversion factor dt p / dt have this qualitative behavior. In other words, some of the speed-up 
is accounted for in the parametric equations and the rest is contained in dt p /dt. 



3.2. Comparison with Physical Orbit Solutions 



One way to compare the dynamics of the physical problem and the spirographic ap- 
proximation is to write the effective equations of motion in the same form. The equations 
of motion for the orbit problem itself are given in the text above. The equation of motion 
for the spirograph equivalent system can be written in the form 



d0_ = 1 + (1 - 2(3/a)e 



1/2 



(28) 
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where £ = (x 2 + y 2 ) 1 ^ 2 is the radial coordinate in the orbital plane. Here we consider the 
system to be spirograph equivalent if the orbits have the same turning points (£1,^2) and 
the same turning angle (A0, which sets (3 /a). 

Next we note that the equation of motion for both physical orbits and the spirographic 
orbits can be written in the form 



d9 _ 1 



(£-6) (6-0 1/2 9(0, (29) 



where the functions g(£) are slowly varying over the radial range of the orbits £1 < £ < £2- 
We denote the functions g(£) as distortion functions because they determine the manner in 
which the orbits are distorted from an elliptical shape. For physical orbits in a Hernquist 
potential and the corresponding spirographic approximation, the distortion functions g(£) 
take the form 

rq l + 1/2 ftfr + (1 - 28/a)e ,, n , 

W(f)= l-^J and M0= __ F . (30) 

Notice that for Keplerian orbits, one obtains the same general form as equation (29) with 
<?(£) = \foT e = V£i£2 — constant. The function g phys reduces to this form in the limit a — > 1. 
The corresponding function g spi approaches the Keplerian form in the limit of low energy 
(e — > 0) and high angular momentum (q — > q max ). 

Figures 5 and 6 show how closely the actual (physical) orbits can be approximated 
by the spirographic treatment developed here. Figure 5 shows the radial excursion of an 
orbit, with the dimensionless radius plotted as a function of angle, for a physical orbit and 
the corresponding spirographic orbit. The two curves are nearly indistinguishable, which 
vindicates the approximation for this system (e = 0.90, q/q ma x = 0.75). 

Next we need to quantify the difference between the physical orbits and the spirographic 
approximation to the orbits over the entire parameter space. We can characterize the avail- 
able parameter space in terms of the quantities (e,q/q max ), where both variables have the 
range < c,q/q max < 1. The difference between the two orbits at a particular radius £ can 
be measured by the quantity 

A 9 _ flphys(Q ~ &pi(0 / 31 v 
9 ~ #phys(0 

where the distortion functions g are defined above. The difference between two orbits can be 
measured by the root mean square (RMS) of the difference (Ag)/g averaged over an entire 
orbit, i.e., over the range of radii £1 < £ < £2- For a benchmark comparison, we find the 
portion of parameter space for which the RMS difference is less than 10%. This value of 
10% is arbitrary. However, in many applications, such as the LMC orbit considered in the 
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following section, the observational errors are typically in the range 10 - 20%. The result 
is shown in Figure 6. For the region in the plane above the solid curve, the orbits can be 
represented as epicycloids (spirograph patterns) with an effective error less than 10%. Figure 
6 also shows the portion of the e — q plane for which the orbits are close to Keplerian in the 
same sense. The Kepler approximation is also characterized by a function g (see above), and 
the dashed curve marks the locus of points for which the Keplerian approximation differs 
from the physical orbit by 10%. The Keplerian approximation is thus valid only for the 
portion of parameter space to the left and above the dashed curve. The bottom portion 
of the e — q plane shown in Figure 6 corresponds to low angular momentum and hence 
nearly radial orbits. The portion of parameter space below the dot-dashed curve can be 
modeled, again within 10%, as radial orbits (for which we have analytic solutions - see §3.4). 
Notice that for sufficiently deep orbits (e>0.9), the combination of the spirographic and 
radial approximations is sufficient to describe the orbits; this portion of parameter space 
corresponds to < 0.073. Finally, the location of the orbit of the LMC is shown in the 
diagram, and is found to be squarely within the spirographic regime (as discussed in §4). 

Another way to assess the validity of the spirographic approximation is to see how well 
the orbits conserve energy and angular momentum. Conservation of energy is built into the 
approximation (through the conversion factor dt p /dt) and is satisfied exactly. However, the 
spirographic orbits follow slightly different paths (than the physical orbits) and experience 
small variations in angular momentum over the course of an orbital cycle. Within the 
spirographic formalism, the angular momentum in dimensionless form is given by j spi = 
\xv y — yv x \dt p /dt, which can be written 



Jspi 



1 + e 



— e 



i/ 2 7 2 (a/jg _ i) _ ( a _ + 7 ( a - (3)( a /(3 - 2) C os{ai p /P) 
[ 7 2 (a//3 - l) 2 + (a - f3) 2 - 2 7 (a - P)(a/(3 - 1) cos(at p /(3)] 1/2 ' 



If the spirographic approximation were exact, this expression would be constant. In general, 
the approximation scheme leads to a small variation (Aj)/j over the course of an orbit. 
The dotted curve in Figure 6 shows the trajectory in parameter space for which the angular 
momentum has an RMS deviation of 10% over the course of an orbit. In the region above 
the curve - most of parameter space - the variation in angular momentum is less than 10%. 
Notice that the requirement of angular momentum conservation, as enforced here, is less 
stringent than that of orbital shape, as measured by the distortion functions (compare the 
dotted and solid curves in Fig. 6). Since conservation of angular momentum is not precisely 
satisfied, the spirographic formalism does not provide an exact solution to any approximate 
physical problem; on the other hand, it does provide an extremely good approximation to 
the exact physical problem. 
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3.3. Orbital Decay under a Prictional Force 

Within the spirographic approximation, we can calculate how orbits change under the 
action of a frictional force. As an example, we consider the frictional force f per unit mass 
to be proportional to the velocity so that 

f=-^v, (33) 

where the leading coefficient is constant and T has units of time (this form applies to the low 
speed limit of dynamical friction, e.g., see Binney & Tremaine 1987). The torque exerted 
on the orbiting body is given by r x f and the time rate of change of the specific angular 
momentum j in terms of spirographic elements becomes 

| = -"f {VI 2 ~ (a - P) 2 + (a- PM V - 1) cos[(l + v)Q}^ , (34) 

where rj = (a — 0)f (3 and where the conversion factor di p /dt is given by equation (27), 
including the fiducial time scale to- The change in angular momentum over each half orbit 
can be written 

Aj = dt^ t = - r f J*'" dt p { V1 2 -(a- P) 2 + (a- /J) 7 fa - 1) cos[(l + V )t p ]} , (35) 

where we have changed the variable of integration (to parametric time) in the second equality. 
This procedure allows us to evaluate the integral to obtain 

Aj = - r f[ m 2 -(a-f3) 2 }((3n/a). (36) 

For comparison, in the limit of circular orbits (where a, (3 — > and r] — > 1) the change in 
angular momentum reduces to Aj = —{r- 2 s /T)^ 2 {n/2). This latter form is that expected 
from applying a constant torque over a specified time (a half orbital period). The difference 
between the circular expression and equation (36) shows how the geometry of the orbit affects 
the decay of angular momentum. Within the spirographic approximation, this expression 
for the change in angular momentum is exact. 

Next we consider the loss of orbital energy. The work done on the orbiting body by the 
frictional force over a half orbit leads to a loss of energy, which can be written 

/r 2 r n /2 rjf 
f • vdt = -f j o dt { 7 V + (a - P) 2 - 2 1V (a - (3) cos[(l + v)i P ]}(^) 2 ■ (37) 

The integral can be simplified by changing the integration variable to parametric time. In 
this case, however, one factor of di p /dt is left over. Invoking the Mean Value Theorem, we 
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can pull the extra factor out of the integral and the expression for AE becomes 

If we could identify the correct value of (di p / dt) , then the above expression would be exact 
within the spirographic formalism. In practice, we have to settle for an approximation; the 
mean value {dt p /dt) ~ (/3vr/ari/ 2 ) will suffice for most applications. In summary, the loss of 
energy per half orbit becomes 

AE » -=^-[tY + (a - /3) 2 ](/Wa) 2 • (39) 

J- Ti/2 

In addition to the approximation made in using the spirographic form (which typically leads 
to errors of a few percent), this expression for energy loss contains a few percent error 
incurred in the evaluation of the integral. 



4. THE ORBIT OF THE LARGE MAGELLANIC CLOUD 

As an application of the solutions found here, we consider the orbit of the LMC in 
the potential of the Galaxy. Using observational data in conjunction with the treatment 
developed above, we find the orbital elements of the LMC motion. The closest equivalent 
spirographic approximation reproduces the physical orbit to a precision of ~ 6 — 7%, and we 
thereby obtain a parametric (and hence analytic) description of the orbit. 

For this demonstration, we take the following observed quantities as given: The observed 
velocities of the motion of the LMC in galactic coordinates (van der Marel 2002) are as 
follows: radial velocity v r = 84 ± 7 km/s, tangential velocity vq = 281 ± 41 km/s, and 
total velocity vt = 293 ± 39 km/s. The quoted errors lie in the range 8 — 15% and result 
in uncertainties of comparable order in the derived quantities of the orbit. The distance 
to the LMC is now reasonably well known, with r\ mc ps 50 kpc (e.g., van der Marel 2002). 
We also need to specify the total mass M in contained within the LMC orbit; using velocity 
data from a large number of Milky Way satellites, including the LMC, Kochanek (1996) 
finds M in = 5 ± 1 x lO n M . These quantities are sufficient to specify the orbit within the 
framework developed here, i.e., assuming that the potential of the Galaxy can be modeled 
as a Hernquist potential over the range of radii probed by the LMC orbit. For completeness, 
we note that in order to reproduce the proper Galactic rotation curve at the solar circle, 
we would need to include additional disk/bulge components which are neglected here (and 
which introduce another source of uncertainty in the orbital parameters). 
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Conservation of energy implies 

\v 2 T = * [-e + 7— -7 — ] ~ (207km/s) 2 , (40) 

* 1 + Slmc 

where £i mc = r\ mc /r s . The scale ^0 of the potential can be written in terms of known 
quantities, 

, ta= GM^ = „ (207km/s)2 . (41) 

'"s '"lmc Slmc Slmc 

We can solve the above two equations for the dimensionless energy e to obtain 

_ 1 + £lmc — 6mc//^ 1 
' " (1 + 6mc) 2 ^ (l+6mc) 2 ' 



(42) 



where we have defined = 2GM in /(ri mc t>f,). The current observational estimates suggest 
that /i ~ 1, which leads to the second approximate equality. Using the definition of q, we 
can write 

^(2GM in /r lmc )- 1 l - 6m / , 2 » — 7- 6m / , 2 , (43) 



9 = 



(l+6mc) 2 /X (1 + 6 



lmc ) 



where we have used the observational results quoted above to obtain the numerical coefficient. 

If the scale length r s of the Galaxy and mass M in are known, equations (42) and (43) 
specify the energy and angular momentum of the LMC orbit. The working estimate for M in 
uses data from many satellites of the Milky Way (Kochanek 1996) and is reasonably secure 
(so that /x ~ 1 ± 0.2). We estimate r s as follows: Numerical simulations of the formation 
of dark matter halos show that the density distributions approach a nearly universal form 
(NFW) when properly scaled. A recent numerical study (Navarro et al. 2004) indicates that 
the halo density profiles should be scaled so that they line up at r_ 2 , the radius at which the 
logarithmic slope of the density profile s = d\ogp/d\ogr = —2. The paper also shows that 
for halos with galactic masses (rotation velocities ~ 200 km/s), the radius r_ 2 = 18 — 28/2T 1 
kpc. For the Hernquist density profile used here, the logarithmic slope s = —2 occurs at 
£ = 1/2 so that r s = 2r_ 2 (note that this relation is different for the original NFW profile 
where r_ 2 = r s ). Using h = 0.7, we estimate the scale length of the Galaxy to be r s m 65 kpc. 
As a result, £i mc = r\ mc /r s 0.77 and our estimates for the energy and angular momentum 
of the LMC orbit are 

e ~ 0.32 and q ^ 0.134 (q/q max w 0.69) . (44) 

The solid hexagon in Figure 6 marks the location of the LMC orbit in the plane of possible 
orbits. Notice that the orbit falls within the region where the spirographic approximation is 
valid. The observed quantities used to estimate the orbital parameters are uncertain at the 
~ 10% level, so the resulting orbital elements are subject to comparable uncertainties. 
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Using the values of e and q found here, we can specify the other parameters of the LMC 
orbit. The inner turning point r\ ~ 46 kpc, the outer turning point r 2 ~ 114 kpc, and 
the effective semi-major axis of the orbit r* = r s £* ~ 82 kpc. With these estimates for the 
turning points, the generalized eccentricities of the LMC orbit are e = (£ 2 — £i)/(2£*) ~ 0.41 
ande= (^2 — ^1) / (^2+^1) ~ 0.43 The turning angle of the half orbit is given by A6/71 « 0.706, 
so the orbit is far from closing. The total radial period r 1.75 Gyr, which is comparable 
to the period of the equivalent circular orbit r c = 27r(r lmc 3 /G'M in ) 1 / 2 ks 1.5 Gyr. 

For the estimates of e and q found here, the RMS deviation of the physical orbit from 
its spirographic approximation is ~ 6.6%. Because this difference is much smaller than the 
uncertainties in the orbit due to observational error and systematics, there is no practical 
difference between the actual physical orbit and its spirographic approximation. For com- 
parison, the RMS deviation of the physical orbit from the Keplerian approximation is much 
larger (about 17%). The spirographic orbital elements (in physical units) are thus a = 48 
kpc, /3 — 14 kpc, and 7 = 80 kpc. The corresponding spirographic orbital eccentricity e s 
= (a — /3)/7 ~ 0.425. Figure 7 shows an overlay of the physical orbit and its approximate 
form given by the spirographic formalism. The curves are shown in the orbital plane and 
the phases of the orbits are aligned. The close agreement between the two orbits argues that 
the spirographic approximation is adequate for modeling the LMC orbit. 

In the discussion thus far, we have assumed that the scale length r s and mass M in 
of the Galaxy within the radius of the LMC are known. The possible variation of these 
quantities provides an important source of uncertainty in the orbital elements. For example, 
if we assume that the scale length r s is completely unknown, then equations (42) and (43) 
constrain the LMC parameters to lie along the curve defined by 

q « 0.92(1 -v^) 3 e~ 1/2 . (45) 

If we allow the scale length r s to vary over the range 50 < r s < 80, bracketing the estimated 
value, the allowed orbits in the e — q plane lie on the locus of points marked by open squares 
in Figure 6. Similarly, we can vary the mass of the Galaxy within the current LMC position 
by varying the mass parameter ji in the above equations (/i = M in /(5 x 1O 11 M ). If we let 
the enclosed mass vary over its allowed range /i — 1 ± 0.2 (Kochanek 1996), the orbit lies 
along the locus of points marked by open triangles in Figure 6. 

This analysis implies a lower limit on the mass M in . In order for the LMC to be in a 
bound orbit, e > and equation (42) implies that 

f 1 > T, 7 7 = • I 46 ) 

(i + 6mc) n mc + r s 

The scale length r s cannot be arbitrarily large. If we adopt a conservative upper limit on r s 
based on the results of numerical studies (e.g., Navarro et al. 2004, Bullock et al. 2001), then 
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r s < 100 kpc and hence /i > 1/3. The corresponding mass limit becomes M in > 1.7 x lO n M . 
This limit is consistent with previous estimates of Kochanek (1996) and others (van der Marel 
2002; Wilkinson k Evans 1999; Lin, Jones, & Klemola 1995). 

Next we apply the frictional formulas developed in §3.3 to the orbit of the LMC (see 
also Murai & Fujimoto 1980). The above analysis determines the values for the spirographic 
parameters a, f3, and 7. Using these results in equation (36), we find that the change in 
angular momentum over a half orbit is given by Aj « — 3.12r;?/T. The angular momentum 
itself is given by j = r s y/2qWo. Putting these results together with the parameters of the 
LMC orbit, we find that the fractional change in angular momentum per orbit is given by 



The time scale T can be evaluated using previously known results (Binney & Tremaine 1987). 
The new result found here is that the change in angular momentum (per orbit) obtained 
using the true orbital shape (Aj oc [rj^ 2 — (a — f3) 2 ]f3 /a) is larger than that of the equivalent 
circular orbit (Aj oc 7 2 /2) by 31%. 

5. ORBITAL PROPERTIES FOR GENERAL MASS DISTRIBUTIONS 

The discussion thus far has focused on orbits in the Hernquist potential and spirographic 
approximations to those orbits. In this section we show that a wider class of potentials 
leads to orbits with similar shape and can be adequately modeled using the spirographic 
approximation. Toward this end, we prove a series of results that apply to all physical 
potentials (§5.1 and 5.3) and find the region of parameter space for which the spirographic 
approximation is valid for a collection of particular potentials (§5.2). 



This section presents a set of general constraints on orbit shapes for any extended mass 
distribution. This set of constraints not only delimits the allowed orbital paths, but also 
shows that spirographic curves (epicycloids) provide reasonable approximations to the orbits 
for general extended mass distributions. The potential for any such mass distribution allows 
for two and only two turning points for bound orbits (Contopoulos 1954). This constraint 
limits the orbital path to lie within the annulus £1 < £ < £2- We also know that the turning 
angle for a half orbit is confined to the range n/2 < A9 < n (Contopoulos 1954). To 
proceed further, we assume that the turning points and the turning angle (£1, £ 2) and A9) 




(47) 



5.1. General Constraints on Orbital Shape 
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are given, and find constraints on the path taken as the orbit travels from (£ 2 , 0) to (£1, A9) 
in the orbital plane. To illustrate these constraints, which are general, we plot one particular 
example in Figure 8. Here, we use turning points £i = 0.284, £ 2 = 0.893, and a turning angle 
A9 = 2.01 (for a Hernquist profile, these values correspond to the choices of dimensionless 
energy e = 0.5 and angular momentum q/q max = 0.5, but the constraints derived below hold 
for any physical potential). 

One can show that the orbits contain no inflection points (S. Tremaine, private com- 
munication). A related constraint can be found by using the definition of the turning angle 
and finding an upper bound on the magnitude of the derivative \d£/d9\. Since the potential 
is monotonic, ip < = e + g/£i 2 over the radii of interest, and this bound takes the form 

<k =([(7d 2 -l] 1/2 . (48) 



d9 

The integral of this quantity produces a curve of minimum radius £ for a given angle 9, and 
this curve has the form 

£(^) = £ 1 csc[^ + sin- 1 (£ 1 /6)]. (49) 

This function starts at the outer turning point and reaches the inner turning point "as soon 
as possible", in the sense that any real orbit must travel inward more slowly as a function 
of angle 9. This argument produces a complementary limiting form - the orbit could stay 
at the outer turning point for as long as possible and then plunge inward (along a path of 
maximum \d£/d9\) just in time to reach the inner turning point at 9 — A9. The true orbit 
is confined to lie between the two aforementioned curves, which are shown as the solid lines 
in the illustrative example of Figure 8. 

We can also consider the slowest possible turning of the orbit. The turning angle A9 < ir 
(Contopoulos 1954) and only attains the value of 7r in the Keplerian limit. The proof of this 
result shows that the magnitude of the derivative \d£/d9\ is always greater than that of 
the Keplerian orbit. As a result, the orbit must fall (from the outer turning point towards 
the inner turning point) faster than the Keplerian orbit with the same turning points, and 
the true orbit must lie within the Keplerian ellipse. As before, this argument provides a 
complementary constraint: If we construct a Keplerian orbit with the same turning points 
that reaches the inner turning point at 9 — A9, and trace the orbit backwards out to the 
outer turning point, the physical orbit must lie outside the curve. The physical orbit must 
thus lie between these two ellipses, which are shown as dashed curves in Figure 8. 

In general, as the density profile become less centrally concentrated, the turning angle 
decreases. In particular, A9 = n in the limit of a point mass (the Keplerian limit), and the 
turning angle A9 = n/2 for a uniform density distribution (p = constant). Since physical 
orbits are confined to have turning angle A9 > ir/2, we argue that the uniform density orbit 
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falls toward the inner turning point faster than any other case. The uniform density model 
thus provides a benchmark for comparison and its orbits are determined by the function 
/ c (0 defined by 

/c(0 = (£ 2 -6 2 )(6 2 -£ 2 )- (so) 

We want to show that physical orbits are bounded by the orbit of the uniform density model, 
with the turning points £1,^2 chosen to coincide with those of the physical orbit. Consider 
the difference in angle between an arbitrary physical orbit [determined by /(£)] and the 
constant density benchmark case [determined by / c (01- We define the difference in turning 
angle at a given radius to be = /) — f c ) so that 

Now we look for the maximum angular difference. The derivative dQ/d^ = if and only if 
/(£)/? = fc(0/(£i&) 2 a ^ some intermediate point £x < £ < £ 2 (recall that f/q = /cACiG) 2 
at the turning points £1,^2 by construction). Next one can show that /(£)/<? = f c (Q / (£,i£,2) 2 
at an intermediate point only if the second derivative of the physical potential d 2 ip/d^ 2 is less 
than that of the constant density potential (since d 2 ijj c /d^ 2 < for the constant density case, 
"less than" means larger in magnitude and negative). Finally, we argue that the constant 
density potential has the greatest negative curvature of any physically relevant potential, 
so that the derivative dQ/d£, 7^ 0. This result implies that the difference in turning angle 
between any physical orbit and the orbit of the constant density potential must be monotonic. 
By inspection we find that Q is monotonically increasing. As a result, physical orbits are 
confined to lie outside the limiting curve found by integrating the orbit of a constant density 
potential (see eq. [50]), i.e., 

0{t) — — — -sm — — 7; 7T- . (52) 

As before, a complementary bound exists. Since the orbit cannot turn faster than equation 
(52), one can construct a solution in which the orbit stays at the outer turning point as long 
as possible and then follows a function of the form (52) to the inner turning point. Physical 
orbits must lie within such a curve. This argument leads to the pair of dotted curves shown 
in Figure 8. 

The above constraints limit the trajectory taken by an orbit in any extended mass 
distribution. Each of the constraints leads to a pair of bounding curves, with the true 
physical orbit confined to lie between each pair of curves. The resulting combination of 
constraints implies an orbit much like that of the epicycloid (spirographic) curve. This set 
of constraints is illustrated in Figure 8 for one particular case, but this class of constraints 
holds for all orbits. 
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5.2. Allowed Parameter Space for Specific Potentials 

Another way to illustrate the usefulness of the spirographic approximation is to ex- 
plicitly calculate the allowed regions of parameter space for a collection of extended mass 
distributions. In addition to the Hernquist profile discussed above, we calculate orbits for 
three additional potentials and find the regimes of parameter space for which the spirographic 
approximation is valid. The first additional model is the NFW profile, where the density 
distribution and potential can be written in the dimensionless form 

P= and ^ = ^ln(l+0- (53) 

Next we consider a density distribution that is more concentrated and has the form 

"=WUT7¥ and '^(TW (54) 

We denote this case as the 3/2 model (from the power-law slope of its density profile in the 
inner limit). Like the Hernquist model, these potentials reach a finite central value and we 
have scaled the dimensionless fields so that ^(0) — 1. As a result, the energy e is confined to 
the range < e < 1. For a given energy, each potential allows a maximum value q max of the 
angular momentum parameter (q max corresponds to a circular orbit) and the q parameter is 
confined to the range < q/q ma x < 1- 

We also consider the Jaffe model, where the density profile has the form of a singular 
isothermal sphere for small radii, but falls more quickly at large radii and reaches finite total 
mass (Jaffe 1983). These profiles have the form 

P = ^T^ ^d ^ = V ln(l±i). (55) 

Unlike the previous cases, this potential does not reach a finite value at the center and hence 
the dimensionless energy e does not have a finite range. The other potentials described 
above are defined over the energy range < e < 1 and allow orbits over the full radial extent 
< £ < oo. In order to compare the Jaffe model with these other potentials, we adopt the 
value Vo = 1/5; with this normalization, the orbits for e = 1 have an effective semi-major 
axis pa 0.0041 (instead of — > 0). 

For each potential defined above, we have calculated the orbits. For each physical 
orbit, we obtain the turning points £ 2 ) an d turning angle A9, and use the results of §3 to 
define a spirographic approximation to the orbital path. The departure of the true orbit from 
its spirographic approximation is then calculated using the method of distortion functions 
developed in §3.2. For these potentials, Figure 9 shows the region of parameter space for 
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which the spirographic approximation to the orbits is valid, with an accuracy of (at least) 
10%. The allowed region of the e — q plane lies above the various curves: Hernquist profile 
(solid curve), NFW (dashed curve), the 3/2 profile (dotted curve), and the Jaffe model (dot- 
dashed curve; keep in mind that the parameter space for the Jaffe model extends beyond the 
portion of the e — q plane shown here). At low energies corresponding to large radius orbits, 
the spirographic approximation has a limited range of validity for all of the potentials. At 
higher energy, the regime of parameter space for spirographic orbits depends on the degree 
of central concentration of the density profile. The Hernquist and NFW profiles have the 
same density dependence for small £ (p ~ £ _1 ) and similar ranges of validity. In general, 
the region of validity of the spirographic approximation shrinks as the power-law index of 
the inner density profile increases (compare the previous models with the 3/2 model and the 
Jaffe model). This finding, that less centrally concentrated mass distributions give rise to 
orbits that are more spirographic in shape, can be understood from analytical considerations: 
In the limit of a uniform density distribution, the turning angle of the orbits has the value 
A6 = tt/2 and the corresponding distortion function is identical to that of the spirographic 
approximation (see eq. [30] in the limit f3 — » a/2). For completness, we note that this 
analysis only considers the shape of the orbits. In order to provide a complete description of 
the dynamics, one must include the transformation between parametric time and physical 
time (the analog of eq. [27]), and this transformation depends on the potential. 

Each potential also has a region of parameter space for which the Keplerian approxi- 
mation is valid (these regions are not delimited here). The density profiles for the Hernquist 
and Jaffe models approach the form p ~ £~ 4 at large radii and hence reach the Keplerian 
regime for sufficiently large £ (low energy e). In contrast, the enclosed mass for the NFW 
profile diverges logarithmically and hence the Keplerian approximation is never applicable. 
The 3/2 model lies in between these two cases, with a narrow regime of parameter space (at 
the low e edge of the plane) for which the Keplerian approximation is valid. 

5.3. Definitions of Orbital Eccentricity 

The definition of orbital eccentricity necessarily contains an ambiguity for extended mass 
distributions. Specifically, we show here that the two alternate definitions of eccentricity 
e = (C,2 — £i)/2£* and e — (£2 — £i)/(£2 + £1) are equivalent if and only if the potential is 
Keplerian. The first definition e measures eccentricity relative to the circular orbit given 
by £* and has a maximum value e max < 1- The second definition e is not centered on the 
circular orbit, but attains the full range < e < 1. The first half of the argument is well 
known - the two definitions of eccentricity are the same for a Keplerian potential. 
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To complete the argument, we suppose that the claim is true so that e = e. Then 
£1 + £2 = 2£*, where is the orbital radius of the circular orbit (for a given e), and £ x and 
£2 are the turning points. For a general potential ip{£,), the turning points are given by the 
zeroes of the function /(£), defined by 

/(0 = -^ 2 + eV(0-9- (56) 

The radius of the circular orbit is given by the condition df /d£ = 0. In order for £1 + £2 = 
2£* to be valid for all orbits, the function /(£) must be symmetric with respect to the point 
£ = £*, i.e., + rr) = — x) for all x G [0, £*]. If the function /(£) is symmetric, 
all derivatives of odd order must vanish at £ = £*. The first derivative df/d£ vanishes by 
definition. The third derivative must also vanish and is given by 

-de = ^W + ^W + Q ^ = - (57) 

This condition does not depend on e or q, but only on the form of the potential. Furthermore, 
the possible range of e and q allow for all values of to be realized. In order for the 
third derivative of /(£) to vanish for all orbits, equation (57) must hold for all values of £. 
This differential equation is second order in the function U = dip/d^ and has two linearly 
independent solutions: U a = 1 / 1; 2 &nd Ub = l/£ 3 - The first solution corresponds to ipA ~ l/£ 
which is the Keplerian case. The second solution leads to ip B ~ l/£ 2 , which implies a negative 
mass density, and is thus unphysical. The original differential equation (57) is third order 
and has a third solution, ipc = constant, which is not of interest. Thus, we have shown that 
if e = e, then the potential ip — l/£ and is Keplerian. In order words, the validity of the 
condition £1 + £2 = 2£* requires a Keplerian potential. 



6. CONCLUSION 

This paper has explored orbital trajectories for the Hernquist potential and other ex- 
tended mass distributions. Our results can be summarized as follows: 

[1] The first set of results quantitatively determines the orbits for the Hernquist po- 
tential. The analysis is straightforward. The orbits have a rosette shape (Fig. 1), as was 
well-known previously in qualitative terms (Binney & Tremaine 1987). The resulting orbits 
can be described by their turning angles A6* and half periods Ti/ 2 for a given value of energy 
e and angular momentum q (see Figs. 2 and 3). In analogy to the Kepler problem, we have 
defined osculating orbital elements (Table 1). The Hernquist profile allows for a number 
of results to be derived analytically. We have found the energy and angular momentum as 
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function of the turning points (eq. [6]), the maximum angular momentum for a given energy 
(eq. [9]), a radial scale that plays the role of the semi-major axis (eq. [10]), the analog 
of Kepler's law (§2.2), an analytic form for the turning angle in the limit of circular orbits 
(eq. [16]), and an expression for the radial period in the limit of low angular momentum (eq. 
[17]). This latter expression provides a good approximation to the radial period for all orbits 
(Fig. 3). We have also constructed a fitting formula (eq. [11]) that specifies the turning 
angle A6>(e, q) as a function of energy and angular momentum over most of parameter space. 

[2] This paper shows that the physical orbits of the Hernquist profile can be modeled 
to a good approximation by epicycloid curves, more commonly known as spirograph pat- 
terns. These curves have analytic solutions, which can be written in the parametric form of 
equation (20). These parametric equations allow for a completely analytic description of the 
orbits, and we have developed this approximation in some detail. The transformation be- 
tween spirographic orbital elements and physical parameters is given by equations (21-22). 
The effective eccentricity of a spirographic curve is given by equation (23). Most impor- 
tantly, the spirographic approximation faithfully reproduces the shape of the physical orbits 
to an accuracy better than 10% over most of parameter space (Fig. 6). The spirographic 
approximation conserves energy exactly (using the transformation of eq. [27] to relate phys- 
ical time to parametric time) and conserves angular momentum to a greater accuracy than 
it reproduces the orbital shape. 

[3] As a demonstration of the efficacy of this approach, we have found the shape of 
the orbit of the Large Magellanic Cloud as it traces through the halo of our Galaxy. If 
we model the Galactic potential with a Hernquist potential, the dimensionless energy and 
angular momentum of the LMC orbit are (e, q/q m ax) = (0.32, 0.69) and the orbit takes the 
form shown in Figure 7. The spirographic approximation conserves angular momentum to 
an accuracy of 1% and reproduces the shape of the physical orbit to an accuracy of 6 - 
7%. These levels of error are much smaller than both the observational uncertainties in 
the problem and the errors incurred in approximating the Galaxy as a smooth, spherical, 
extended mass distribution. As a result, the analytic (spirographic) approximation provides 
a good working model for the orbit (with essentially no loss of accuracy). 

[4] In addition to considering orbits in the Hernquist potential, this paper argues that 
the results found here - in particular the spirographic approximation to the orbits - apply to 
general extended mass distributions. Toward this end, we have presented general constraints 
on the orbital shape (§5). Figure 8 shows that the orbits are restricted to be relatively close 
to the shape of an epicycloid for any extended mass distribution. Figure 9 shows that the 
spirographic approximation is valid over much of the e — q plane for a collection of specific 
potentials (including the NFW profile). We have also shown that the definition of orbital 
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eccentricity necessarily contains an ambiguity: For a given mass distribution, the eccentricity 
can either be defined to be symmetric with respect to a radial scale that plays the role of 
the semi-major axis, or the it can be defined solely in terms of the turning points (see eq. 
[12]); the two definitions coincide only for a Keplerian potential (§5.3). 

Many astrophysical systems involve orbits in extended mass distributions and an ana- 
lytic approximation to the orbits is often a useful tool. This paper develops the spirographic 
approximation, which can be used to model orbits in a wide variety of contexts. One such 
example, that of the LMC orbit, has been studied in this paper, but many additional ap- 
plications of this formalism remain. Possible future topics include orbits of galaxies within 
their clusters, orbits of dark matter particles for dark matter detection strategies, globular 
cluster orbits in galactic potentials, orbits of small halos as they merge into larger halos 
during structure formation, orbits of young stars in their birth aggregates, and many others. 

In addition its direct applications, the approach developed here can be generalized in 
a number of ways. This work has focused on the Hernquist potential, although some com- 
parison with other models has been included (see Fig. 9). Other cases should be considered 
and compared with the spirographic approximation, including the HH model in which two 
different Hernquist potentials are nested together (Ciotti 1996), models with varying degrees 
of central concentration (Tremaine et al. 1994), and the inclusion of a point mass at the 
origin (to model the central black holes found in many galactic centers). For power-law 
potentials, the turning angles have already been computed (Touma & Tremaine 1997) and 
the spirographic approximation can be readily implemented. The orbits in this paper are 
confined to a single plane, for spherical potentials, but the spirographic approximation can 
be implemented for two dimensional orbits in nonspherical potentials. 

We note that the Hernquist potential and alternate forms (e.g., the NFW profile or 
the Jaffe model) are just models for more complicated physical structures. In the case of 
galactic halos, for example, the true potential must have corrections due to dumpiness and 
hence small scale inhomogeneities; a galactic disk (or multiple disk components) and hence a 
large scale quadrapole moment; perturbations from nearby large galaxies (e.g., Andromeda); 
and other complications. Many of these effects enter at the ~10% level and imply 10% 
corrections to the orbital shape calculated from a spherical potential. This paper shows 
that over most of the relevant parameter space, the spirographic approximation leads to 
orbital shapes that agree with those of the Hernquist profile (and others - see Fig. 9) to an 
accuracy of better than 10%. Thus, in practice, the spirographic approach can be used as 
an approximation to the orbits without losing accuracy (compared to numerically calculated 
orbits) while providing an analytic description of the dynamics. 
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Fig. 1. — An example of an orbit in the Hernquist potential with energy e = 0.10 and angular 
momentum parameter q = q ma xl^- The orbit does not close and can be characterized by the 
angle through which the orbit turns during one radial oscillation. One part of the orbit is 
highlighted (with the open square symbols). 
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Fig. 2. — The turning angles for one half of an orbit. The quantity A9 is the angle through 
which the orbit turns as the radius varies from the inner turning point to the outer turning 
point. The turning angle A0 for a half orbit is plotted here as a function of q/q m ax, where 
Qmax is the maximum angular momentum accessible for a bound orbit at a given energy. 
The various curves correspond to dimensionless energy values e = 0.01, 0.1, 0.3, 0.5, 0.7, and 
0.9 (from top to bottom in the figure). Notice that A6 — > n/2 in the limit q — » 0, but, in 
general, A9 ^ it in the limit of circular orbits (q — > q max ). 
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Fig. 3. — The range of orbital half period. The quantity T\/2 is the time required for the 
orbit to travel from the inner turning point to the outer turning point. Here, the half period 
is normalized by two functions, first the factor e 3 / 2 , which represents the variation expected 
in the Keplerian limit, and by the factor [cos -1 yfi — y/e(l — e) 1 / 2 ], which gives the correct 
orbit time in the limit of zero angular momentum (q — > 0). The half period is plotted as a 
function of q/q max , where q max is the maximum angular momentum accessible for a bound 
orbit at a given energy. The various curves correspond to dimensionless energy values e = 
0.01, 0.1, 0.3, 0.5, 0.7, and 0.9 (from bottom to top in the figure). 




Fig. 4. — Allowed range of parameter space for spirographic orbits (epicycloids) to approxi- 
mate physical orbits. The drawing point parameter 7, which plays the role of the semi-major 
axis, has been set equal to unity. The parameter a represents the radius of the larger circle 
and j3 is the radius of the smaller circle. The value of (3 must lie below the dashed line 
(where [3 = a/2) so that the turning angle of a half orbit is greater than ix/2. The (3 value 
must also lie above the dotted line (where f3 = a — 7) so that the inner turning point is 
positive. Physical orbits are confined to the region below the dashed line and above the 
dotted line. The inset diagrams show representative orbital shapes for the delimited regions 
of parameter space. The upper left portion of the diagram, above the solid line where a = (3, 
is not allowed {(3 < a by definition). The upper right portion of the diagram (above the 
dashed line and below the dotted line) allows for an alternate representation of the physical 
orbits (see text). 
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Fig. 5. — Comparison of the radial coordinate as a function of turning angle for the physical 
orbit problem in a Hernquist potential (solid curve) and the equivalent spirograph orbit 
(dashed curve). The two systems are chosen to have the same radial turning points (£i and 
£2) and the same angular displacement per half orbit A9 1.84 (Ad/n 0.585). These 
values correspond to the energy e = 0.90 and the angular momentum variable q/q m ax = 
0.75 for the Hernquist potential. According to the method of distortion functions developed 
in the text, the two orbit shapes differ by 2.2%; the spirographic orbit conserves angular 
momentum to an accuracy of 0.3%. For comparison, the shape of the physical orbit differs 
from that of a Keplerian orbit with the same turning points by 
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Fig. 6. — Regions of parameter space for which the Keplerian and spirographic approxi- 
mations are valid for the Hernquist profile. The region of the e — q plane above the solid 
curve corresponds to orbital parameters for which the RMS deviation in the shape of the 
physical orbits from that of the spirographic approximation to the orbits is less than 10%. 
Similarly, the region above the dashed curve is where the RMS deviation of the physical or- 
bits from Keplerian orbits is less than 10%. The region above the dotted curve is where the 
spirographic approximation conserves the angular momentum of the orbit to an accuracy of 
10%. Finally, the region below the dot-dashed curve is where the orbits are radial to within 
10%. The large hexagon at (e, q/q max ) ~ (0.32, 0.69) marks our best estimate for the orbital 
parameters of the LMC. The accompanying (short) curves show the variation of the orbital 
parameters with variations in the Galactic scale radius r s (open squares) and enclosed mass 
M in (open triangles). 
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Fig. 7. — The predicted shape of the LMC orbit in the potential of the Galaxy. The 
coordinates x and y represent the plane of the orbit and are expressed in kpc. The solid 
curve shows the orbital shape for a Hernquist potential using our best estimate for the 
dimensionless orbital energy and angular momentum (e,q/q ma x) = (0.32,0.69). The dashed 
curves shows the spirographic (analytic) approximation to the orbit. The shape of the orbits 
agree to within about 7%, as measured by the distortion functions (eqs. [30, 31]). As another 
measure, the analytic (spirographic) orbit conserves angular momentum to an accuracy of 
1%. For comparison, the observational determination of the orbital parameters is subject to 
10 - 20% uncertainties. 
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Fig. 8. — General constraints on orbital shape for an arbitrary potential. The bold-faced 
curve shows the spirographic orbit for given turning points (marked by the open symbols, 
with radii marked by the inner and outer circles) and a given turning angle for a half orbit. 
The two solid lines represent an upper limit on how fast the orbit can turn, so that the orbit 
must trace a path between the two lines. Similarly, the dashed curves depict a limit on how 
slowly the orbit can turn; this limit corresponds to elliptical orbits and physical orbits must 
fall between these curves as well. The dotted curves represent the constraint that the orbit 
cannot turn any faster than the orbit of the uniform density model (which has A0 = 7r/2 
and all orbits must have A9 > n/2). Physical orbits (for given £i,£ 2) A0) must lie between 
the three pairs of constraints, thereby limiting the orbital shape to be reasonably close to 
that of the spirographic orbit. 
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Fig. 9. — Regions of parameter space for which the spirographic approximation is valid for a 
collection of potentials. The region of the e — q plane above the curves corresponds to orbital 
parameters for which the shape of the physical orbits deviates from that of the spirographic 
approximation by less than 10%. The potentials considered here are the Hernquist potential 
(solid curve), the NFW profile (dashed curve), the 3/2 model with density profile p ~ 
£-3/2^ _|_ ^/£)- 4 (dotted curve), and the Jaffe model (dot-dashed curve). The Jaffe model 
differs from the other cases in that the potential does not reach a finite central value and 
hence the plane depicted here does not represent all of its parameter space. 



